setwd("C:/Users/86139/Desktop/R作业")
library(Seurat)
library(dplyr)
#数据读取####
pbmc.data <- Read10X(data.dir = "filtered_gene_bc_matrices/hg19")
pbmc <- CreateSeuratObject(counts = pbmc.data, 
                           project = "pbmc3k", 
                           min.cells = 3, 
                           min.features = 200)
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc)
#特征gene选择####
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
#scale
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
#PCA降维
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
#细胞聚类####
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)
pbmc <- RunUMAP(pbmc, dims = 1:10)
LabelClusters(DimPlot(pbmc, reduction = "umap"),id = 'ident')
#细胞注释####
new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "B", "CD8 T", "FCGR3A+ Mono", 
                     "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
#差异分析####
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
top2_markers <- pbmc.markers %>%
               group_by(cluster) %>%
               top_n(n = 2, wt = avg_log2FC)
print(top2_markers$gene)
DotPlot(pbmc, features = top2_markers$gene) + RotatedAxis()
FeaturePlot(pbmc,features = c("CCR7","LEF1","FOLR3","S100A12","AQP3","CD40LG","LINC00926"), pt.size = 0.5) + NoLegend()
VlnPlot(pbmc, features = c("CCR7","LEF1","FOLR3","S100A12","AQP3","CD40LG","LINC00926")) + NoLegend()
